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Abstract. In the context of the Instantaneous Normal Mode approach, the spectrum 
of the Hessian of Hamiltonian is a key quantity to describe liquids behaviour. The 
determination of the spectrum represents a major task for theoretical studies, and 
has been addressed recently in various works. In this work a new approach for the 
analytic computation of the Hessian spectrum is presented. The one dimensional case 
for a system of particles interacting via a purely repulsive potential at low density is 
Q ■ analyzed in details and the spectrum is computed exactly also in the localized sector. 

Finally, the possible extensions of the method are discussed, together with a comparison 
with different approaches to the problem. 
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1. Introduction 

In this work we describe an approach we have recently developed to compute analytically 
the so called Instantaneous Normal Mode (INM) spectrum of a liquid system |l] . That 
is, more precisely, the density of eigenvalues of the Hessian of the Hamiltonian, averaged 
y\f ' over the equilibrium distribution. 

The general frame where our computation acquires particular meaning, and where 
a great number of recent analysis have been performed, is the Instantaneous Normal 
Modes (INM) approach ||. The main idea of this approach is that liquids are 'solid- 
like' at short times and that liquids' dynamics thus correspond to vibrations about some 
equilibrium positions with periodic jumps into new local minima ||. 
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In this context, a crucial quantity is the typical spectrum of the Hessian of the 
Hamiltonian, which describes the structure of the energy landscape around the typical 
configurations. This is the INM spectrum and the related eigenvectors are the so called 
Instantaneous Normal Modes. 

Recently there have been many attempts to relate quantitatively the knowledge of 
the INM spectrum with measurable quantities, such as the velocity-velocity correlation 
function, more complicated self-correlation functions, and even the diffusion properties 
||. The main idea is to generalize the standard harmonic analysis used for solids, taking 
into account the finite hopping rate from one local minimum to the other j|, |J. For 
many systems simulations have shown that the predictions given by the INM analysis 
are in very good agreement with the results obtained directly from molecular dynamics. 
Simulations also show that the INM spectrum of a liquid system always displays both 
positive and negative eigenvalues (i.e. real and imaginary frequencies), not only in 
the liquid, but also in the supercooled phase as the glass transition is approached. 
While the interpretation of the positive modes seems straightforward (they represent 
harmonic vibrations in well defined wells), this is not the case for the negative modes. 
The contributions to these modes come from regions of negative curvature in the phase 
space sampled by the system at equilibrium. Indeed the system does not remain forever 
in a local minimum of the energy surface, but jumps into new regions over certain time 
scales, in this way exploring an enormously complicated landscape. 

Unfortunately, it is not possible to distinguish in a simple way the modes related 
to real barriers (which can be used to give an estimate of the hopping rates, see ref. ||) 
from the ones related to anharmonic deformations of the local landscape. Recently it has 
been proposed || that a crucial difference exists between localized and extended negative 
modes: the localized modes involve a finite number of particles and can be associated to 
local barriers hopping; the extended modes involve an extensive number of particles and 
can be associated to structural rearrangements. Thus these two kinds of modes should 
be related to different diffusional processes. Simulations seem to confirm this view for 
some fragile models: in || it is shown that the number of extended negative modes 
goes to zero as a particular temperature is approached and that this temperature is the 
same as the temperature T c where the Mode Coupling Theory (MCT) would predict 



a dynamical transition ||1Q| . The temperature T c represents an important reference 
value for fragile glasses: above T c MCT successfully predicts the observed dynamical 
behaviour of the system, and the viscosity has a non-Arrhenius temperature dependence; 
below T c MCT breaks down and the viscosity behaviour becomes Arrhenius-like. It is 
commonly believed that T c coincides with the crossover temperature postulated long 
ago by Goldstein [[J, below which hopping processes become relevant. The result of || 
shows that, for the particular system studied there, T c can be found via the spectrum 
analysis, looking for the temperature where the extended negative modes disappear. It is 
not clear how general this phenomenon is, but it anyway indicates that the localization 
properties of the INM represent an useful tool, and can be used to give at least an 
estimate of the temperature T c . 
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For all these reasons, the analytic determination of the INM spectrum represents 
a vital task for any theoretical study of liquids. Indeed there have been in the last 
years various attempts to perform such a computation and some important steps in this 
direction have been done in M and [|7]] . In the following we will try to outline what the 
standard procedure to compute analytically the INM spectrum is, and we will briefly 
mention what are the main assumptions and strategies adopted in these previous works 
|| [5J. This will help us in introducing our new approach, commenting on what are 
the differences in perspective and procedure, and what are the objectives we propose to 
address. 

2. General procedure 

Let us consider a system of N interacting particles with Hamiltonian 

N 

H = Y,V(r kl ), (1) 

k>l 

where V(r) is a two-body potential. The Hessian matrix A is defined by, A^ = d^d\H , 
with fi,u = 1, . . . , d, being d the dimension of the space. The general form of A with 
respect to particle indices is 

JV 

Akl = —Jki + &kl 2^ Jki , (2) 

i 

where J kl = J{r kl ) = V"(r kl )r H r H + ^p(l - hihi), hi being the versor along the 
inter-particle k, I distance. The diagonal term of A is a consequence of the translational 
invariance of the system, which requires J2k ^-ki — 0. 

The standard procedure to compute the INM spectrum is to consider the well 
known relation between the density of eigenvalues of the matrix A and the correspondent 
resolvent operator G(A) = (Al — A) -1 : 

D A (X) = lim -— Im Tr G(A - it) . (3) 

This formula provides the spectrum of a single matrix A, that is of the Hessian matrix 
evaluated in a particular particles configuration. To obtain the INM spectrum -D(A), 
one has to average Eq. (^ over the Boltzmann equilibrium distribution. In this way 
the computation of the INM spectrum is reduced to the computation of the average 
diagonal element of the resolvent matrix (Gu(\ — ie)). The standard procedure is 
at this point to use an integral representation for the resolvent: {Ga(X — it)) = 
j? f d(f>i . . . d(pN <t>1 exp[— |0(A — A + ie)<f)] . In this way some new variables, the internal 
fields 4>i, enter in the computation, beside the space coordinates r^ included in the explicit 
expression of the Hessian A. A further step consists in adopting a replica trick to bring 
the normalization factor Z at the numerator: -| = lim n ^ Z n ~ x . In this way we finally 
get: 

(G«(A -ie)) = lim fd$i (<Pl) 2 Q($i) (4) 



A new approach for the analytic computation of the Instantaneous Normal Modes spectrum 4 
where now the fa are vectors in a n-dimensional space. The function Q is given by 

fl(S) = fdfa... d$ N (e-^ x ~ A+i <>*> , (5) 



where the integration is performed over all, but the i, internal fields. 

Of course the difficult task is to compute the function Q. In [§,0], more or less 
explicitly, the authors assumed a Gaussian shape for £l(<j>), and then computed self- 
consistently its variance using a sort of generalized liquid theory (we will come back 
to this point later). Thus, they assumed a quite rough approximation for the general 
4> dependence of Q, but took into account accurately the many particles contribution 
terms included in the definition ([5]). Our approach has been precisely the opposite 
one. As you will shortly see, we did not assume any a priori form for Q(4>), but we 
disregarded many particles contribution terms. As a result, we have an approach which 
is much simpler to deal with, but still allows for non trivial spectral properties. In || [7]], 
excellent results were obtained for the density of eigenvalues, but, on the other hand, 
the computational procedure was too complicated to reasonably look at the localization 
properties of the eigenf unctions. Besides, a Gaussian shape for Q is anyway too simple 
to detect non trivial localization properties. On the contrary, our simpler approach is 
more suitable for an eventual study of localization. 

3. Our approach 

We consider a liquid system at low density, that is precisely in the physical context 
where we can expect many particle correlations not to be important. At low density 
the most part of the particles will be very far one from the other and we can imagine 
that, picking out at random a couple of particles % and j, the interaction Vfaj) between 
them will be very low, practically zero. This will be true also for the Hessian element 
Aij between the two particles, since the Hessian is much more short ranged than the 
potential itself. The whole matrix A will therefore have the most part of its elements 
equal to zero and very few ones different from zero: it will be, in other terms, a diluted 
matrix. In this context, the main idea underlying our approach is to model the matrix A 
as a random diluted matrix and to use all the techniques developed for random matrix 
ensembles to compute the INM spectrum. 

As a first step in this program, we must find out what is the probability distribution 
of the matrix A, or, which is the same, of the matrix J. This distribution is 
naturally induced by the equilibrium probability over the positions of the particles 
via expression (0). Consistently with our low density approximation, we can assume 
that the probability distribution P[J] is factorized into the individual probabilities of 
the particles pairs. In this way the elements of J (but not of A) are independently 
distributed, i.e. 

N 

p[j) = Hp(j kl ). (6) 

k>l 
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It is clear that with (|6]) we are disregarding three-particles correlations. Thus, (|6|) 
becomes a reasonable assumption when three-particles correlations are not important, 
for example at low densities. 

Once assumed this f actor ized form for P[J], we can express the pair-probability 
p(Jki) as, 



p(Jki) = J dn . . . dr N e pH 5(Jki - Jim)) 



= ^Jdr k dn g {2 \r kl )5(J kl - J(r kl )) (7) 

where g^(r) is the two-particles correlation function and p is the average density. At 
this point there are various ways to exploit equation ([<]). One can for example insert 
in (|7|) the numerical values obtained for g^ 2 \r) by numerical simulations, or by some 
classical liquid theory approaches. Alternatively, and this is the simplest possibility, we 
can adopt a low density expansion for the two points correlation function. In this way, 
g^ 2 \r) = exp(— (3V(r)), at the first order of the virial expansion. This approximation 
is consistent with our previous assumption @ where we disregarded three-particles 
correlations, and, as we shall see, it enables to perform completely analytically the 
calculations. 

In the general frame we have described, we finally have the Hessian matrix A, and 
its distribution, as given by (|7|). This is enough to start a random matrix computation. 
However, the algebra is still very complicated: the matrices involved (A and J) are 
actually tensors with particle and space coordinates indices. Before dealing with this 
general case, we decided to test our method under the simplest possible conditions, that 
is when the particles live in one dimension. This case is conceptually analogous to the 
three dimensional one, with the advantage of a simpler algebra. 

4. The one dimensional case 

In one dimension the explicit expression of the matrix J is much simpler, since only the 
longitudinal part of @ survives, and we have J(r) = V"{r). To obtain the distribution 
p(Jki) we have to solve equation (|7|). In the light of the greatest simplicity we consider 
a soft-spheres potential V(r) = l/r m . Thus, from eq.(|7]) we have, 

1 e~P Jb 1 
P(J) -Wl^^N q{J) • (8) 

with (3 = (3 [m(m + l)] m + 2 , b = m/(m + 2) and c = l/(m + 2). (From now on 
we indicate with J an individual element of the matrix J). For realistic values of m 
(typically m = 12) the parameter b is very close to one. Therefore, we will directly set 
b = 1 in p( J) in order to simplify our calculation. We will show in the discussion of our 
results that the actual spectrum is very weakly dependent on this approximation. 

As it stands the distribution p(J) is not normalizable, but we can regularize it in 
the following way. Let us put an IR cut-off f, by setting V(r) = for r > f, and let 
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tj = V(f). We obtain in this way a regularized form of the pair probability: 

p v (J) = 6(J) + i (q(J) 0(J - V ) ~ 8{J) j™ df q{ J')) , (9) 

where q(J) is defined in equation (|[). A few comments about equation ([|): 

• the distribution p v (J) is diluted as we expected, since the probability of finding an 
element of the matrix J equal to zero is of order one, while the probability of finding 
one element larger than 77 is of order 1/N. 

• Pri(J) explicitly depends on the value r\ of the cut-off. However, if we consider 
a generic function f(J) and compute its average over p v (J), we get: (f(J)) = 
f,^° dJ q(J) If (J) — /(0)]. Thus, if the function / is differentiable in 0, its average 
value has a well defined limit when rj — > 0. This means that in our computation, 
after averaging over the Hessian distribution, we can safely take the limit r\ — > 0, 
recovering the original problem without cut-off. 

At this point we finally have a well defined random matrix problem: we have 
an ensemble of matrices A, and their distribution, as given by equation (^). We can 
therefore try to apply random matrix techniques to compute the density of eigenvalues 
D(X). 

Before going on with the computation, we would like to add some more general 
remarks on the one dimensional case. 

First of all, we note that the low density approximation we made when disregarding 
three-particles correlations in ([]) is, at fixed density and temperature, the less 
appropriate the lower the dimension. For d — 1 the geometrical constraints on the 
particles are much stronger and indeed the true Hessian matrix has a band structure. 
For this reason our analysis has to be regarded more as a training example for the three 
dimensional case rather than a predictive computation for a real one dimensional system 
(our one dimensional treatment is very similar to what, for an electronic band structure 
problem, would be an s-band computation |TT|). 

Secondly, in the specific example we are dealing with we have chosen a soft-sphere 
potential V(r). It is easy to see that in this case the Hessian A is a positive defined 
matrix which therefore has a positive defined density of eigenvalues. This is not the case 
for a three dimensional system, where the role of the negative modes is an important 
issue. However, the qualitative shape of the spectrum is actually very similar to what 
found, for example, in simulations on three dimensional systems [|13 |. Moreover, the 



eigenf unctions exhibit non trivial localization properties in the tails, thus providing a 
good context where testing a localization analysis procedure. Finally, to avoid confusion 
in the future, we note once again that the matrix J is not a nearest-neighbors matrix, 
but a diluted one. This is why, even in one dimension, we find both localized and 
extended states (and not only localized ones |fT2|). 



Let us now proceed with the computation. The procedure is very similar to the 
one outlined in the previous section for the standard computation of the INM spectrum: 
one has to relate D(\) to the resolvent operator G(A), and the diagonal element of the 
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resolvent to the one-particle function Q. The main difference is that now all the averages 
involved in the definition (|5]) of fl are not averages over the Boltzmann distribution, but 
over the distribution p v (J). 

The crucial point is, of course, the computation of fl(<pi). This can be done in 



different ways, generalizing some random matrix computations [Q to the particular 
case of the distribution p v ( J) |I[. In our case it is convenient to write fl as (the particle 
index i is, from now on, understood) fl((p) = e _ 2 A< ^ +9(<f>) anc [ i 00 k f or a se if consistent 
equation for the exponent g(<p). We note that g(<fr) measures how much the function fl 
is Gaussian: a non-quadratic shape of g implies a non-Gaussian fl. 

We will not enter into the details of how obtaining the self consistent equation (the 
interested reader is referred to |l|), but simply give the result: 

g {$) = f <0 e -W' 2 +s(<?) f d j g (j) L-iJ(<?-<?') 2 _ t l (10 ) 



This equation has still to be averaged over the distribution q(J). This can be done 
exactly with a few algebraic tricks, and finally we get 

f°° A 

g(x) = g(e ln/4 x) = K x (x) - x / dy K 2 (x, y) exp(i-?/ 2 + g(y)) (11) 

Jo I 

where x = |0|, and Ki(x) and K 2 (x,y) are expressed in an analytic form |l[. 

First of all, it is possible to check analytically that asymptotically g(x) ~ x 2c , 
thus proving that fl is definitely not a Gaussian function. Besides, we have been able 
to numerically solve the equation for g(x) without any further approximation. Indeed, 
eq.flTTD has the form of a fixed-point equation and can be solved numerically by iteration, 
discretizing the function g and the kernel K on a lattice. 

Once obtained g for a given value of A, it is possible to compute the spectrum, 
using backward the relations between fl and the resolvent, and finally between the 
resolvent and D(\) (see previous sections). The results are shown in Fig. la, where we 
have plotted the INM spectrum D as a function of A, for m = 12. The spectrum has 
positive support and it depends on the scaled inverse temperature (3 in the expected 
way: for low temperatures (high p) the collisions among particles are weaker, so that 
the spectrum is peaked on lower values of the eigenvalues. On the other hand, the tail 
for large A is larger at higher temperature. The behaviour in the right tail is of the form 
D{X) ~ e~ aX , where a is an increasing function of f3. 

As previously stated, this one dimensional case represents for us, first of all, a way 
to test our analytic procedure. A crucial task is therefore to check whether the result we 
have found is correct. To this aim we have done extensive numerical simulations. Once 
drawn a matrix J with probability @, we build A and diagonalize it numerically. Since 
the spectrum has huge tails for large eigenvalues, it is convenient, in order to compare 
simulations with analytic results, to consider the probability distribution n of /i = In A, 
that is vr(/i) = D(e^)e^. In Fig. lb we plot 7r(/i) as obtained from the analytic form of 
D(X), together with the one obtained from numerical simulations. The two curves are in 
excellent agreement confirming the validity of our result. Besides, we show in the inset 
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of Fig. lb the numerical spectrum obtained with the original value of b = m/(m + 2). 
The result justifies the sensibility of the approximation b ~ 1. 

In the introduction we mentioned that an important issue in the context of the 
INM calculations is the analysis of the localization properties of the negative modes. As 
we have previously said, even if in this one dimensional case the spectrum is positive 
defined, the tails exhibit non trivial localization properties. From a numerical point of 
view an important quantity to investigate localization properties is the average inverse 
participation ratio Y(X a ) = YhLi ( w a) > where a = 1 ... N is the eigenvalue index and 
w a = [( a \i)} 2 i s the weight of site i in the eigenfunction \X a ). In Fig. 2 we plot Y as a 
function of A, as obtained via numerical diagonalization. It is clear from the figure that 
there are two localizations edges, separating a central region of extended eigenvalues, 
from the tails where localized states are present. Note that for A — ^ the inverse 
participation ratio goes to one and this corresponds to a single particle which happens 
to be nearly decoupled from the rest of the system. On the other hand, the localized 
states of the right tail correspond to pairs of very strongly interacting particles and this 
naturally leads to a inverse participation ratio equal to 1/2. 

Our aim for the future is to find analytically the two localization edges revealed by 
the numerical analysis. Our idea is to generalize to our case the methods used for the 



Bethe Lattice [|I(J and we will address this problem in a future work |[20|| . 

5. Comparison with other approaches 

In order to better understand what are the differences between our approach and the 
ones of [0, [7] , it is convenient to go back once again to the general procedure outlined in 
section II, and, in particular, to expressions (§) and (H). A crucial observation made in 
is that the function Q(<p) is proportional to the one-particle correlation function of a 
generalized liquid system, whose coordinates are X t = (r*, (pi) and whose Hamiltonian 
is given by: 

- PH GEN = ]T U(Xi) + 2 W{X h X,) , (12) 

i i<j 

where U(Xi) = -|A (fa ■ fy and W(X t ,Xj) = -/3V(ry) + Mj {4>i • <?j)- Indeed it 
is simple to check that £l(<f>) = (^1^5(0 — (pi)) gen = s (A)/p, where s(X) is the 
one-particle correlation function of the generalized liquid system, and p is the average 
density of the original system. Given this analogy, it is possible to compute Q using the 
standard liquid theory |17| applied to this generalized liquid system. This is precisely 



what the authors of did: they assumed a Gaussian shape for Q, and then computed 
the variance using a renormalized Mean Spherical Approximation |17|, [18], |19| . 

An alternative route is the following. We can use a Hypernetted Chain (HNC) 



approximation [|I7], [18], M to obtain the one-particle correlation function of the 
generalized liquid. The HNC equations can be deduced from a variational principle 
where the free energy of the system is written as a functional of the one-particle and 
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two-particles correlation functions 
form: 



I8|. The variational free energy has the following 



[3F = ]-J dXdX's(X)s(X') £ (2) (X, X') log(£ (2) (X, X')) + \- J dXdX' 



s{X)s{X') l-g^(X,X f )+(3W(X,X') g^(X,X' 



(2), 



dXs(X)U(X) 



+ J dXs(X)[\og(s(X)) 



^ ~ iv Tr 



log(l + sh) — sh+ -(hshs) 



(13) 



where s(X) and g^'(X,X') are, respectively, the one-particle and the two-particles 
correlation functions, h(X, X') = g( 2 >(X,X') — 1 and in the trace term all the products 
are convolutions. Self-consistent equations for s(X) and g^ 2 '(X, X') are obtained 
variationally from (0). These equations are in general not easy to solve, but one can 
start facing them perturbatively. For example, we can look at a low density situation 
where, in a first approximation, the trace term will not contribute. In this case, 
the variational equation for the two-particles correlation function gives immediately 
g( 2 '(X, X') = e~ /3W ^ x,x '^ , which is nothing else than the first order of the virial 
expansion. The equation for the one-particle correlation function is, on the other hand, 
less trivial. In terms of fi(0) = s(X)/p it gives 



log(fi(0)) + -A0 2 = l<lo'9.(,>')i> 



dre -PV{r) e \J{r){$-$>? _ N 



(14) 



The interesting fact about equation fll~4]) is that it is exactly the same self-consistent 
equation that we get in our low density approach presented in the previous section. To 



see this, we note that p f dre ^^ej 



J(r)(</>-</>') 2 



fdJq(J) 



e-2 



JiWr _ 1 



+ iV. If we 



use now the previously adopted notation Q((fi) = exp(— 1/2 X(p 2 + g((p)) and we insert 
it in (|TJ[) we finally get precisely equation ( [T0| ) for g(4>). 

This result shows that our simple low density approach can be seen as the first 
order approximation of a more complicated generalized HNC approach. We understand 
now in a deeper way what is the origin of the approximations we made and also what is 
the route we have to follow to improve our calculation. Indeed, we see that, if we want 
to go beyond the low density approximation, we have to consider further terms in the 
trace appearing in ( |i~3l) and compute consequently the self-consistent equations ]20 . 



6. Conclusion 



We have presented a new approach for the analytic computation of the INM spectrum. 
At present we have successfully applied this approach to the case where particles live in 
one dimension, and at low density. This simple case represents the ideal context where 
understanding the potentiality of our method and outlines the procedure one has to 
follow in the more general cases. 

We expect this approach to give good results in the realistic three-dimensional 
situations, where the geometrical constraints are much less important than in one 
dimension. A strong indication in this direction is provided in Fig.2b, where we plot the 
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spectrum obtained from a Montecarlo simulation for a soft-sphere mixture. The thick 
curve is obtained considering, for each sampled configuration, the real Hessian matrix 
A of the system. The points, on the other hand, are obtained considering a 'scrambled' 
Hessian matrix where, for each sampled configuration, a new J matrix is built by mixing 
at random the elements Jij of the real J. This scrambling procedure destroys the three- 
particles correlations and is therefore equivalent to our approximation (|J). The figure 
shows clearly that the scrambled approximation works very well for T = p/3 1 ^ 4 = 0.2. A 
more systematic analysis |20| with Montecarlo simulations indicates that good results 



are obtained also at much higher values of T suggesting that our analytic approach, 
even at the simplest first step where three-particles correlations are discarded, should 
describe quantitatively well the real spectral properties of a three-dimensional system 
on a wide range of temperatures and densities. 
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n 0.20 




Figure 1. (a) The INM spectrum Dasa function of A for different values of the 
scaled temperature (3, m = 12 and e = 0. The plot is in log-linear scale. In the inset it 
is shown D(\) for (3 = 1.5 in linear-linear scale. The spectrum vanishes at A = 0. (b) 
Numerical simulations vs. analytic solution. We plot here the probability distribution 
7r(/i), with /i — In A. Both the curves correspond to a probability distribution with 
b = 1. N = 600 and 77 = 10~ 4 . (3 = 1 and m = 12. In the inset, on the same 
scale, we compare the analytic result for b = 1 with the simulations performed with 
b = m/(m, + 2). 
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Figure 2. (a) Inverse participation ratio Y as a function of the eigenvalue A at 
different values of N. r/ = 10~ 4 , (3=1 and m = 12. The thick curve is the spectrum 
D(A) at the same values of the parameters, (b) Montecarlo simulation for a three- 
dimensional soft-sphere mixture, at T = 0.2. The thick curve represents the real INM 
spectrum. The points correspond to the spectrum obtained using a scrambled Hessian 
matrix. 



